options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)
install.packages('cmdstanr',repos=c('https://mc-stan.org/r-packages/',getOption('repos')))
library(cmdstanr)
check_cmdstan_toolchain(fix=T)
install_cmdstan(cores=2)
cmdstan_version()
restart R
options(mc.cores=parallel::detectCores())
make stan file '---.stan'
Bernoulli distribution
data{
int N;
array[N] int y;
}
parameters{
real<lower=0,upper=1> p;
}
model{
p~beta(1,1);
y~bernoulli(p);
}
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}
data=list(N=9,y=sample(0:1,9,replace=T))
mdl=cmdstan_model('./ex0-0.stan')
#mdl$exe_file()
mle=mdl$optimize(data=data)
## Initial log joint probability = -7.2058
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 5 -6.18265 0.000131667 2.68779e-08 1 1 8
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.3 seconds.
mle
## variable estimate
## lp__ -6.18
## p 0.56
mcmc=mdl$sample(
data=data,
#seed=1,
#chains=4,
iter_sampling=500,
iter_warmup=100,
#thin=1,
refresh=1000
)
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
# or
#mcmc=goMCMC(mdl,data)
mcmc
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -8.07 -7.80 0.72 0.30 -9.52 -7.58 1.00 1029 1211
## p 0.55 0.55 0.14 0.14 0.30 0.77 1.00 748 1291
#or
#mcmc$summary()
drw=mcmc$draws('p')
par(mfrow=c(4,1),mar=c(1,5,1,1))
drw[,1,] |> plot(type='l',xlab='',ylab='p')
drw[,2,] |> plot(type='l',xlab='',ylab='p')
drw[,3,] |> plot(type='l',xlab='',ylab='p')
drw[,4,] |> plot(type='l',xlab='',ylab='p')
par(mar=c(3,5,3,3))
par(mfrow=c(1,2))
drw %>% hist()
drw %>% density() %>% plot()
seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(4,1),mar=c(1,5,1,1))
drw[,1,] |> plot(type='l',xlab='',ylab=pr)
drw[,2,] |> plot(type='l',xlab='',ylab=pr)
drw[,3,] |> plot(type='l',xlab='',ylab=pr)
drw[,4,] |> plot(type='l',xlab='',ylab=pr)
par(mar=c(3,5,3,3))
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}
mcmc$metadata()$stan_variables
## [1] "lp__" "p"
mcmc$metadata()$model_params
## [1] "lp__" "p"
seeMCMC(mcmc)
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -8.07 -7.80 0.72 0.30 -9.52 -7.58 1.00 1029 1211
## p 0.55 0.55 0.14 0.14 0.30 0.77 1.00 748 1291
saveRDS(mdl,'bin_mdl.rds')
saveRDS(mcmc,'bin_mcmc.rds')
mdl=readRDS('bin_mdl.rds')
mcmc=readRDS('bin_mcmc.rds')
library(qs) #faster
qsave(mcmc,'bin_mcmc.qs')
mcmc=qread('bin_mcmc.qs')
n=30
tb=tibble(x=runif(n,0,9),y=rnorm(n,x,1))
tb
## # A tibble: 30 × 2
## x y
## <dbl> <dbl>
## 1 2.83 2.45
## 2 4.26 3.76
## 3 7.60 8.19
## 4 6.92 7.07
## 5 5.66 4.49
## 6 1.60 1.88
## 7 4.68 3.63
## 8 8.39 10.3
## 9 0.615 0.233
## 10 7.94 7.47
## # ℹ 20 more rows
f0=formula('y~x')
f0
## y ~ x
lm(f0,tb)
##
## Call:
## lm(formula = f0, data = tb)
##
## Coefficients:
## (Intercept) x
## -0.29 1.04
X=model.matrix(f0,tb) # get explanatory variables from what in formula
X
## (Intercept) x
## 1 1 2.833
## 2 1 4.264
## 3 1 7.599
## 4 1 6.918
## 5 1 5.661
## 6 1 1.600
## 7 1 4.682
## 8 1 8.389
## 9 1 0.615
## 10 1 7.937
## 11 1 7.109
## 12 1 5.520
## 13 1 6.692
## 14 1 8.331
## 15 1 2.101
## 16 1 0.841
## 17 1 0.279
## 18 1 3.446
## 19 1 3.392
## 20 1 0.111
## 21 1 6.372
## 22 1 4.384
## 23 1 8.172
## 24 1 5.698
## 25 1 2.587
## 26 1 3.794
## 27 1 7.949
## 28 1 8.567
## 29 1 7.475
## 30 1 3.166
## attr(,"assign")
## [1] 0 1
tb=tibble(x1=runif(n,0,9),
x2=runif(n,0,9),
c1=sample(c('a','b','c'),n,replace=T),
y=rnorm(n,x1-x2+(c1=='b')*5-(c1=='c')*5,1))
tb
## # A tibble: 30 × 4
## x1 x2 c1 y
## <dbl> <dbl> <chr> <dbl>
## 1 4.34 8.77 b 0.607
## 2 7.94 7.61 c -6.19
## 3 2.28 8.46 c -10.1
## 4 3.25 4.62 a 0.0277
## 5 0.743 1.11 b 5.26
## 6 3.05 6.28 a -2.81
## 7 0.371 0.190 b 4.75
## 8 0.273 0.823 a -1.39
## 9 3.65 7.06 a -3.82
## 10 5.29 0.805 c -0.213
## # ℹ 20 more rows
f0=formula('y~x1')
f0
## y ~ x1
lm(f0,tb)
##
## Call:
## lm(formula = f0, data = tb)
##
## Coefficients:
## (Intercept) x1
## -4.809 0.991
model.matrix(f0,tb)
## (Intercept) x1
## 1 1 4.341
## 2 1 7.937
## 3 1 2.280
## 4 1 3.253
## 5 1 0.743
## 6 1 3.045
## 7 1 0.371
## 8 1 0.273
## 9 1 3.649
## 10 1 5.294
## 11 1 8.480
## 12 1 1.910
## 13 1 2.685
## 14 1 8.004
## 15 1 5.883
## 16 1 3.892
## 17 1 5.148
## 18 1 8.230
## 19 1 1.437
## 20 1 1.155
## 21 1 5.810
## 22 1 7.283
## 23 1 3.715
## 24 1 0.403
## 25 1 3.792
## 26 1 7.194
## 27 1 8.346
## 28 1 1.072
## 29 1 7.658
## 30 1 8.515
## attr(,"assign")
## [1] 0 1
f1=formula('y~x1+x2')
f1
## y ~ x1 + x2
lm(f1,tb)
##
## Call:
## lm(formula = f1, data = tb)
##
## Coefficients:
## (Intercept) x1 x2
## 0.909 1.057 -1.380
model.matrix(f1,tb) # get explanatory variables from what in formula
## (Intercept) x1 x2
## 1 1 4.341 8.775
## 2 1 7.937 7.608
## 3 1 2.280 8.455
## 4 1 3.253 4.624
## 5 1 0.743 1.110
## 6 1 3.045 6.284
## 7 1 0.371 0.190
## 8 1 0.273 0.823
## 9 1 3.649 7.058
## 10 1 5.294 0.805
## 11 1 8.480 6.943
## 12 1 1.910 1.387
## 13 1 2.685 8.230
## 14 1 8.004 2.346
## 15 1 5.883 2.871
## 16 1 3.892 6.618
## 17 1 5.148 5.041
## 18 1 8.230 1.213
## 19 1 1.437 3.285
## 20 1 1.155 5.542
## 21 1 5.810 4.366
## 22 1 7.283 2.415
## 23 1 3.715 4.136
## 24 1 0.403 1.211
## 25 1 3.792 5.762
## 26 1 7.194 5.203
## 27 1 8.346 1.504
## 28 1 1.072 7.879
## 29 1 7.658 6.016
## 30 1 8.515 2.861
## attr(,"assign")
## [1] 0 1 2
f2=formula('y~x1+x2+c1')
f2
## y ~ x1 + x2 + c1
lm(f2,tb)
##
## Call:
## lm(formula = f2, data = tb)
##
## Coefficients:
## (Intercept) x1 x2 c1b c1c
## -0.102 1.088 -0.969 5.046 -5.644
model.matrix(f2,tb) # make categorical variable to dummy variable
## (Intercept) x1 x2 c1b c1c
## 1 1 4.341 8.775 1 0
## 2 1 7.937 7.608 0 1
## 3 1 2.280 8.455 0 1
## 4 1 3.253 4.624 0 0
## 5 1 0.743 1.110 1 0
## 6 1 3.045 6.284 0 0
## 7 1 0.371 0.190 1 0
## 8 1 0.273 0.823 0 0
## 9 1 3.649 7.058 0 0
## 10 1 5.294 0.805 0 1
## 11 1 8.480 6.943 0 1
## 12 1 1.910 1.387 1 0
## 13 1 2.685 8.230 0 0
## 14 1 8.004 2.346 0 0
## 15 1 5.883 2.871 0 1
## 16 1 3.892 6.618 0 1
## 17 1 5.148 5.041 1 0
## 18 1 8.230 1.213 0 1
## 19 1 1.437 3.285 0 1
## 20 1 1.155 5.542 0 1
## 21 1 5.810 4.366 1 0
## 22 1 7.283 2.415 1 0
## 23 1 3.715 4.136 0 0
## 24 1 0.403 1.211 0 1
## 25 1 3.792 5.762 0 0
## 26 1 7.194 5.203 0 0
## 27 1 8.346 1.504 1 0
## 28 1 1.072 7.879 0 1
## 29 1 7.658 6.016 0 1
## 30 1 8.515 2.861 0 0
## attr(,"assign")
## [1] 0 1 2 3 3
## attr(,"contrasts")
## attr(,"contrasts")$c1
## [1] "contr.treatment"
library(brms)
rst=brm(f2,
data=tb,
family=gaussian('identity'), # binomial('logit'),poisson('log')
#seed=1,
#chains=4,
iter=500,
warmup=100,
#thin=1,
prior=c(set_prior('',class='Intercept'),
set_prior('',class='sigma'))
)
rst
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: y ~ x1 + x2 + c1
## Data: tb (Number of observations: 30)
## Draws: 4 chains, each with iter = 500; warmup = 100; thin = 1;
## total post-warmup draws = 1600
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.08 0.48 -1.04 0.86 1.00 894 1070
## x1 1.09 0.06 0.98 1.20 1.00 1946 1096
## x2 -0.97 0.06 -1.10 -0.85 1.00 1596 1348
## c1b 5.02 0.44 4.19 5.90 1.01 545 708
## c1c -5.64 0.38 -6.40 -4.93 1.00 605 612
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.86 0.12 0.66 1.13 1.00 1172 1134
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
stancode(rst)
## // generated with brms 2.21.0
## functions {
## }
## data {
## int<lower=1> N; // total number of observations
## vector[N] Y; // response variable
## int<lower=1> K; // number of population-level effects
## matrix[N, K] X; // population-level design matrix
## int<lower=1> Kc; // number of population-level effects after centering
## int prior_only; // should the likelihood be ignored?
## }
## transformed data {
## matrix[N, Kc] Xc; // centered version of X without an intercept
## vector[Kc] means_X; // column means of X before centering
## for (i in 2:K) {
## means_X[i - 1] = mean(X[, i]);
## Xc[, i - 1] = X[, i] - means_X[i - 1];
## }
## }
## parameters {
## vector[Kc] b; // regression coefficients
## real Intercept; // temporary intercept for centered predictors
## real<lower=0> sigma; // dispersion parameter
## }
## transformed parameters {
## real lprior = 0; // prior contributions to the log posterior
## }
## model {
## // likelihood including constants
## if (!prior_only) {
## target += normal_id_glm_lpdf(Y | Xc, Intercept, b, sigma);
## }
## // priors including constants
## target += lprior;
## }
## generated quantities {
## // actual population-level intercept
## real b_Intercept = Intercept - dot_product(means_X, b);
## }
standata(rst)
## $N
## [1] 30
##
## $Y
## [1] 0.6071 -6.1854 -10.1209 0.0277 5.2596 -2.8055 4.7511 -1.3912
## [9] -3.8195 -0.2133 -2.0143 5.8267 -5.3684 6.5560 -1.4250 -8.8055
## [17] 6.1751 1.6742 -7.2418 -10.0584 8.1638 10.1839 -0.2569 -6.6813
## [25] -0.6644 3.0950 11.4923 -13.1887 -2.9032 5.7330
##
## $K
## [1] 5
##
## $Kc
## [1] 4
##
## $X
## Intercept x1 x2 c1b c1c
## 1 1 4.341 8.775 1 0
## 2 1 7.937 7.608 0 1
## 3 1 2.280 8.455 0 1
## 4 1 3.253 4.624 0 0
## 5 1 0.743 1.110 1 0
## 6 1 3.045 6.284 0 0
## 7 1 0.371 0.190 1 0
## 8 1 0.273 0.823 0 0
## 9 1 3.649 7.058 0 0
## 10 1 5.294 0.805 0 1
## 11 1 8.480 6.943 0 1
## 12 1 1.910 1.387 1 0
## 13 1 2.685 8.230 0 0
## 14 1 8.004 2.346 0 0
## 15 1 5.883 2.871 0 1
## 16 1 3.892 6.618 0 1
## 17 1 5.148 5.041 1 0
## 18 1 8.230 1.213 0 1
## 19 1 1.437 3.285 0 1
## 20 1 1.155 5.542 0 1
## 21 1 5.810 4.366 1 0
## 22 1 7.283 2.415 1 0
## 23 1 3.715 4.136 0 0
## 24 1 0.403 1.211 0 1
## 25 1 3.792 5.762 0 0
## 26 1 7.194 5.203 0 0
## 27 1 8.346 1.504 1 0
## 28 1 1.072 7.879 0 1
## 29 1 7.658 6.016 0 1
## 30 1 8.515 2.861 0 0
## attr(,"assign")
## [1] 0 1 2 3 3
## attr(,"contrasts")
## attr(,"contrasts")$c1
## [1] "contr.treatment"
##
##
## $prior_only
## [1] 0
##
## attr(,"class")
## [1] "standata" "list"
plot(rst) #mcmc trace plot
mcmc=as.mcmc(rst,combine_chains=T) #mcmc sampling
x_new=tibble(x1=runif(3,0,9),x2=runif(3,0,9),c1=c('a','b','c'))
x_new
## # A tibble: 3 × 3
## x1 x2 c1
## <dbl> <dbl> <chr>
## 1 3.98 2.88 a
## 2 0.862 2.32 b
## 3 1.13 8.01 c
# bayes credible interval
fitted(rst,x_new)
## Estimate Est.Error Q2.5 Q97.5
## [1,] 1.45 0.309 0.844 2.07
## [2,] 3.63 0.368 2.910 4.37
## [3,] -12.28 0.371 -13.010 -11.55
marginal_effects(rst,effects='x1') |> plot(points=T)
marginal_effects(rst,effects='x2') |> plot(points=T)
marginal_effects(rst,effects='x1:c1') |> plot(points=T)
marginal_effects(rst,effects='x2:c1') |> plot(points=T)
# bayes predicted interval
predict(rst,x_new)
## Estimate Est.Error Q2.5 Q97.5
## [1,] 1.47 0.927 -0.338 3.27
## [2,] 3.60 0.953 1.641 5.57
## [3,] -12.27 0.944 -14.130 -10.32
marginal_effects(rst,effects='x1',method='predict') |> plot(points=T)
marginal_effects(rst,effects='x2',method='predict') |> plot(points=T)
marginal_effects(rst,effects='x1:c1',method='predict') |> plot(points=T)
marginal_effects(rst,effects='x2:c1',method='predict') |> plot(points=T)
code block
add semicolon ";" at the end of sentence
data{
define objects from R as list
}
parameters{
define estimating parameters from MCMC, can't use int type
}
transformed parameters{
define caluculating parameters from other parameters
}
model{
caluculate log liklihood
define stochastic model, ex. y~normal(m,s) or target+=dist._lpdf(y|m,s)
}
generated quantities{
define caluculating quantities from parameters
}
object type and its difinition of single and array
variable
integer int x, array[N] int x, array[N,M] int x
real real x, array[N] real x, array[N,M] real x
for array[N] int/real x
x[i] x_i
x[i1:i2] array (x_i1...x_i2)
for array[N,M] int/real x
x[i,j] x_ij
x[i] array (x_i1...x_iM)
vector
vector[K] x, array[N] vector[K] x
for vector[K] x
x[i] x_i
x[i1:i2] vector (x_i1...x_i2)
for array[N] vector[K] x
x[i,k] x_ik
x[i] vector[k] x[i]
row_vector[K] x, array[N] row_vector[K] x
simplex[K] x x_i[0,1], Σ x_i=1
unit_vector[K] x Σ x_i^2=1
ordered[K] x x_1 < x_2 ... < x_K
matrix
matrix[J,K] x, array[N] matrix[J,k] x
cov_matrix[K] x symmetric, all eigen>=0
corr_matrix[K] x symmetric, all eigen>=0, x_ij [0,1], x_ii=1
for matrix[J,K] x
x[j,k] x_jk
x[j,] row_vector[K] x_j
x[,k] vector[J] x_k
x[j1:j2,k1:k2] matrix[j2-j1+1,k2-k1+1]
for array[N] matrix[J,K] x
x[i] matrix[J,K] x_i
to fasten a loop for doing matrix[J,K] x
for(k in 1:K){
for(j in 1:J){
do x[j,k]
}
}
parameter constraint
real<lower=0,upper=1> x;
real<lower=l,upper=u> x;
real a[N];
real<lower=min(a),upper=max(a)> x;
vector<lower=0,upper=1>[K] x;
using simplex
categorical distribution
y~Cat(h)
i=1~n, k=1~K, y_i[1,K], h_k[0,1], Σ h_k=1
P[y=k]=h_k
data{
int N;
int K;
array[N] int<lower=1,upper=K> y; // R y[n],y_i[1,K]
}
parameters{
simplex[K] h;
}
model{
for(i in 1:N){
y[i]~categorical(h);
}
}
multinomial distribution
y~multi(n,h)
i=1~n, k=1~K, c_i[1:K], y_k=#(c_i=k), y_k[0,n]
h_k[0,1], Σ h_k=1
P[c=k]=h_k
data{
int K;
array[K] int<lower=0> y; // R y=table(factor(c,levels=1:K))
}
parameters{
simplex[K] h;
}
model{
y~multinomial(h);
}
use vector, matrix instead of array to be fast
using array
data{
int N;
array[N] real x;
array[N] real y;
}
parameters{
real<lower=0> s;
}
transformed parameters{
array[N] real m;
for(i in 1:N){
m[i]= ; // caluculation of x[i]
}
}
model{
for(i in 1:N){
y[i]~normal(m[i],s);
}
using vector
data{
int N;
vector[N] x;
vector[N] y;
}
parameters{
real<lower=0> s;
}
transformed parameters{
vector[N] m;
m= ; //caluculation of x
}
model{
y~normal(m,s);
or
target+=normal_lpdf(y|m,s)
}
mixed model
data{
int N;
int K;
array[N] int<lower=1,upper=K> ID;
vector[N] x;
vector[N] y;
}
parameters{
real a0;
real b0;
vector[K] a;
vector[K] b;
real<lower=0> sa;
real<lower=0> sb;
real<lower=0> s;
}
model{
a~normal(a0,sa);
b~normal(b0,sb);
y~normal(a[ID]+b[ID].*x,s); // .* is multiply element by element
}
K dimension multi normal distribution
data{
int N;
int K;
arrat[N] vector[K] y;
}
parameters{
vector[K] m;
cov_matrix[K] cov;
}
model{
y~multi_normal(m,cov);
}
multi regression
data{
int N;
int K;
matrix[N,K] x; // model.matrix x from R
vector[N] y;
}
parameters{
vector[k] b;
real<lower=0> s;
}
transformed parameters{
vector[k] m;
m=x*b;
}
model{
y~normal(m,s);
}
categorical distribution from raw value
data{
int N;
int K;
array[N] int<lower=1,upper=K> y; // R y[n],yi[1:K]
}
parameters{
simplex[K] h;
}
model{
y~categorical(h);
}
c0=c(1,2,3,4)
h=c(0.4,0.3,0.2,0.1)
y=sample(c0,20,h,replace=T)
table(y) |> prop.table()
## y
## 1 2 3 4
## 0.35 0.15 0.45 0.05
data=list(N=length(y),K=length(c0),y=y)
mdl=cmdstan_model('./ex0-1.stan')
mle=mdl$optimize(data=data)
## Initial log joint probability = -37.0527
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 10 -23.2224 0.000411501 1.05642e-05 1 1 13
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
mle
## variable estimate
## lp__ -23.22
## h[1] 0.35
## h[2] 0.15
## h[3] 0.45
## h[4] 0.05
system.time({
mcmc=goMCMC(mdl,data)
seeMCMC(mcmc,ch=T)
})
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.4 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -31.28 -30.89 1.34 1.03 -34.09 -29.87 1.00 891 972
## h[1] 0.33 0.33 0.09 0.09 0.19 0.48 1.00 3305 1692
## h[2] 0.17 0.16 0.08 0.08 0.06 0.31 1.00 1524 1090
## h[3] 0.42 0.42 0.10 0.10 0.26 0.58 1.00 1376 1474
## h[4] 0.08 0.07 0.06 0.05 0.01 0.20 1.00 819 659
## ユーザ システム 経過
## 1.003 0.259 1.290
categorical distribution from frequency
data{
int K;
array[K] int<lower=0> y; // R y=table(factor(c,levels=1:K))
}
parameters{
simplex[K] h;
}
model{
y~multinomial(h);
}
c0=c(1,2,3,4)
h=c(0.4,0.3,0.2,0.1)
c=sample(c0,20,h,replace=T)
table(c) |> prop.table()
## c
## 1 2 3 4
## 0.35 0.30 0.15 0.20
y=table(factor(c,levels=1:length(c0)))
data=list(K=length(c0),y=y)
mdl=cmdstan_model('./ex0-2.stan')
mle=mdl$optimize(data=data)
## Initial log joint probability = -30.3904
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 8 -26.7017 0.000171529 0.000109974 0.9636 0.9636 11
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
mle
## variable estimate
## lp__ -26.70
## h[1] 0.35
## h[2] 0.30
## h[3] 0.15
## h[4] 0.20
system.time({
mcmc=goMCMC(mdl,data)
seeMCMC(mcmc,ch=T)
})
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -33.98 -33.66 1.25 1.06 -36.36 -32.61 1.00 1142 1543
## h[1] 0.34 0.33 0.10 0.10 0.18 0.50 1.00 3089 1743
## h[2] 0.29 0.28 0.09 0.09 0.15 0.45 1.01 2464 1680
## h[3] 0.17 0.16 0.07 0.07 0.06 0.30 1.00 1258 1422
## h[4] 0.21 0.20 0.08 0.08 0.08 0.36 1.00 1575 1682
## ユーザ システム 経過
## 0.907 0.229 1.192
data{
int N;
vector[N] x;
vector[N] y;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
target+=normal_lpdf(y|b0+b1*x,s);
}
generated quantities{
vector[N] ll;
for(i in 1:N){
ll[i]=normal_lpdf(y[i]|b0+b1*x[i],s);
}
}
n=10
x=runif(n,0,9)
y=rnorm(n,x,1)
data=list(N=n,x=x,y=y)
mdl=cmdstan_model('./ex0-3.stan')
mle=mdl$optimize(data=data)
## Initial log joint probability = -30.4205
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 17 -13.284 0.000308314 0.000853684 0.9301 0.9301 24
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
mle
## variable estimate
## lp__ -13.28
## b0 -0.02
## b1 0.98
## s 0.91
## ll[1] -1.84
## ll[2] -1.04
## ll[3] -1.20
## ll[4] -0.90
## ll[5] -1.06
## ll[6] -1.16
## ll[7] -0.85
## ll[8] -1.53
## ll[9] -1.11
## ll[10] -2.59
system.time({
mcmc=goMCMC(mdl,data)
seeMCMC(mcmc,exc='ll',ch=F)
})
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -15.16 -14.82 1.45 1.24 -17.90 -13.55 1.01 479 536
## b0 -0.07 -0.03 0.79 0.68 -1.35 1.09 1.01 384 312
## b1 0.99 0.98 0.15 0.14 0.76 1.23 1.00 426 390
## s 1.23 1.15 0.39 0.32 0.76 1.98 1.00 752 756
## ll[1] -1.93 -1.79 0.64 0.51 -3.22 -1.16 1.00 1265 1409
## ll[2] -1.34 -1.31 0.35 0.34 -1.97 -0.82 1.00 599 635
## ll[3] -1.47 -1.41 0.43 0.37 -2.23 -0.91 1.00 954 1220
## ll[4] -1.22 -1.20 0.31 0.31 -1.77 -0.76 1.00 615 633
## ll[5] -1.30 -1.28 0.30 0.29 -1.83 -0.85 1.00 1010 982
## ll[6] -1.39 -1.38 0.32 0.31 -1.98 -0.91 1.00 701 865
## ll[7] -1.16 -1.14 0.30 0.31 -1.69 -0.72 1.00 716 760
## ll[8] -1.74 -1.61 0.63 0.52 -2.95 -0.98 1.01 470 395
## ll[9] -1.34 -1.31 0.30 0.30 -1.90 -0.90 1.00 679 676
## ll[10] -2.42 -2.22 0.86 0.72 -4.05 -1.37 1.00 1162 1248
## ユーザ システム 経過
## 0.579 0.236 0.831
ll=mcmc$draws('ll') |>
posterior::as_draws_df() |>
select(contains('ll'))
lppd=sum(log(colMeans(exp(ll))))
pwaic=sum(apply(ll,2,var))
waic=-2*lppd+2*pwaic
data {
int N;
vector[N] x;
array[N] int y;
}
parameters {
real b0;
real b1;
}
model {
target+=poisson_lpmf(y | exp(b0+b1*x));
}
generated quantities {
vector[N] ll;
for (i in 1:N) {
ll[i] = poisson_lpmf(y[i] | exp(b0+b1*x[i]));
}
}
n=10
x=runif(n,-1,1)
y=rpois(n,exp(x))
data=list(N=n,x=x,y=y)
mdl=cmdstan_model('./ex0-4.stan')
mle=mdl$optimize(data=data)
## Initial log joint probability = -52.252
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 10 -10.0548 0.00250835 0.000147357 1 1 12
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
mle
## variable estimate
## lp__ -10.05
## b0 -0.30
## b1 1.25
## ll[1] -1.22
## ll[2] -0.39
## ll[3] -0.42
## ll[4] -1.35
## ll[5] -1.13
## ll[6] -1.20
## ll[7] -0.25
## ll[8] -0.98
## ll[9] -1.49
## ll[10] -1.61
system.time({
mcmc=goMCMC(mdl,data)
seeMCMC(mcmc,exc='ll',ch=F)
})
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -11.04 -10.76 0.99 0.71 -12.92 -10.10 1.00 820 802
## b0 -0.44 -0.40 0.42 0.41 -1.16 0.22 1.00 984 1048
## b1 1.31 1.29 0.61 0.60 0.38 2.29 1.01 692 620
## ll[1] -1.38 -1.29 0.34 0.30 -2.06 -1.01 1.00 757 871
## ll[2] -0.40 -0.35 0.24 0.21 -0.87 -0.11 1.01 734 798
## ll[3] -0.42 -0.37 0.25 0.22 -0.90 -0.12 1.01 744 798
## ll[4] -1.50 -1.40 0.25 0.13 -2.04 -1.31 1.00 1557 1343
## ll[5] -1.26 -1.18 0.26 0.21 -1.78 -1.00 1.00 796 915
## ll[6] -1.27 -1.15 0.31 0.20 -1.94 -1.00 1.00 1749 1450
## ll[7] -0.29 -0.22 0.23 0.17 -0.71 -0.05 1.01 697 661
## ll[8] -0.93 -0.89 0.33 0.31 -1.56 -0.47 1.00 1375 1412
## ll[9] -1.70 -1.59 0.55 0.52 -2.75 -1.03 1.01 710 724
## ll[10] -1.90 -1.71 0.53 0.29 -3.01 -1.50 1.00 1242 1295
## ユーザ システム 経過
## 0.460 0.226 0.734
ll=mcmc$draws('ll') |>
posterior::as_draws_df() |>
select(contains('ll'))
lppd=sum(log(colMeans(exp(ll))))
pwaic=sum(apply(ll,2,var))
waic=-2*lppd+2*pwaic